#ifndef I_LIB_IP_I_FILTER_H
#define I_LIB_IP_I_FILTER_H

#include <cassert>
#include "i_integral.h"

namespace idl
{
    /*Compute the Haar wavelet filter response in x direction.*/
    template <typename T>
    inline T i_haar_x(const T *const *integral, 
        int x, int y, 
        int haar_width, 
        int integral_width, int integral_height)
    {
        int haar_radius = (haar_width >> 1);
        T pos = i_integral_rect(integral, x, y - haar_radius, 
            haar_radius, haar_width, integral_width, integral_height);
        T neg = i_integral_rect(integral, x - haar_radius, y - haar_radius, 
            haar_radius, haar_width, integral_width, integral_height);
        return (pos - neg);
    }

    /*Compute the Haar wavelet filter response in y direction.*/
    template <typename T>
    inline T i_haar_y(const T *const *integral, 
        int x, int y, 
        int haar_width, 
        int integral_width, int integral_height)
    {
        int haar_radius = (haar_width >> 1);
        T pos = i_integral_rect(integral, x - haar_radius, y, 
            haar_width, haar_radius, integral_width, integral_height);
        T neg = i_integral_rect(integral, x - haar_radius, y - haar_radius, 
            haar_width, haar_radius, integral_width, integral_height);
        return (pos - neg);
    }

    /*Compute horizontal derivative 1 0 -1 of a scanline l and store in dx.*/
    inline void i_centered_dx_line(const unsigned char *l, int *dx, int n)
    {
        int i, nm1 = n - 1;
        *dx++ = ((l[0] - l[1]) >> 1);
        for (i = 1; i<nm1; i++)
        {
            *dx++ = ((l[i - 1] - l[i + 1])>>1);
        }
        *dx = ((l[nm1 - 1] - l[nm1]) >> 1);
    }

    /*Compute horizontal derivative 1 0 -1 of an image and store in dx.*/
    inline void i_centered_dx(const unsigned char *image, int *dx, int width, int height)
    {
        int row;
        const unsigned char *pi = image;
        int *pd = dx;
        for (row = 0; row < height; ++row)
        {
            i_centered_dx_line(pi, pd, width);
            pi += width;
            pd += width;
        }
    }

    /*Compute vertical derivative 1 0 -1 between two scanlines lu (up) and ld (down), and store in dy,
    all having n entries*/
    inline void i_centered_dy_line(const unsigned char *lu, 
                                   const unsigned char *ld, 
                                   int *dy,
                                   int n)
    {
        int i;
        for (i = 0; i<n; i++) *dy++ = (((*lu++) - (*ld++))>>1);
    }
  
    /*Compute vertical derivative 1 0 -1 of an image and store in dy.*/
    inline void i_centered_dy(const unsigned char *image, int *dy, int width, int height)
    {
        int row, hm1 = height - 1;
        const unsigned char *pi = image;
        int *pd = dy;
        
        i_centered_dy_line(pi, pi + width, pd, width); 
        pi += width; 
        pd += width;/*first row*/
        for (row = 1; row < hm1; row++)
        {
            i_centered_dy_line(pi-width, pi+width, pd, width);
            pi += width;
            pd += width;
        }
        i_centered_dy_line(pi-width, pi, pd, width);/*last row*/
    }

    /*Run (121)>>2 on n entries of s and store result in n (n must be > 1) entries of d (which could be the same as s).
    Notice that s[-1] and s[n] will Not be accessed, instead, we put special cases at the borders to avoid accessing outside line.*/
    inline void i_binomial_h_121(const unsigned char *s, unsigned char *d, int n)
    {
        int i, count, nm1 = n-1;
        unsigned char a, b, c; /*tmp memorty to store intermidiate results*/
        assert(n > 1);
        a = s[1];
        b = *s++;
        count = (nm1 >> 4);
        for (i = 0; i < count; i++)
        {
            /*add 2 to handle roundoff error*/
            c = *s++; 
            *d++ = (a + (b << 1) + c + 2) >> 2; 
            a = *s++; 
            *d++ = (b + (c << 1) + a + 2) >> 2;
            b = *s++; 
            *d++ = (c + (a << 1) + b + 2) >> 2;
            
            c = *s++; 
            *d++ = (a + (b << 1) + c + 2) >> 2;
            a = *s++; 
            *d++ = (b + (c << 1) + a + 2) >> 2;
            b = *s++; 
            *d++ = (c + (a << 1) + b + 2) >> 2;
            
            c = *s++; 
            *d++ = (a + (b << 1) + c + 2) >> 2;
            a = *s++; 
            *d++ = (b + (c << 1) + a + 2) >> 2;
            b = *s++; 
            *d++ = (c + (a << 1) + b + 2) >> 2;
            
            c = *s++; 
            *d++ = (a + (b << 1) + c + 2) >> 2;
            a = *s++; 
            *d++ = (b + (c << 1) + a + 2) >> 2;
            b = *s++; 
            *d++ = (c + (a << 1) + b + 2) >> 2;
            
            c = *s++; 
            *d++ = (a + (b << 1) + c + 2) >> 2;
            a = *s++; 
            *d++ = (b + (c << 1) + a + 2) >> 2;
            b = *s++; 
            *d++ = (c + (a << 1) + b + 2) >> 2;
            
            c = *s++; 
            *d++ = (a + (b << 1) + c + 2) >> 2;
            
            a = b; 
            b = c;
        }
        count = nm1 & 15; /*handle left columns*/
        for (i = 0; i < count; i++)
        {
            c = *s++; 
            *d++ = (a + (b << 1) + c + 2) >> 2;
            a = b; 
            b = c;
        }
        *d++ = (a + (b << 1) + a + 2) >> 2;
    }
    
    /*Run (121)>>2 to n entries in (s1,s2,s3) and store result in d*/
    inline void i_binomial_v_121(const unsigned char *s1, 
        const unsigned char *s2, 
        const unsigned char *s3, 
        unsigned char *d, int n)
    {
        int i, c;
        c = (n >> 4);
        for (i = 0; i<c; i++)
        {
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
        }
        if (n & 15)
        {
            if (n & 8)
            {
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            }
            if (n & 4)
            {
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            }
            if (n & 2)
            {
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            }
            if (n & 1)
            {
                *d++ = ((*s1++) + ((*s2++) << 1) + (*s3++) + 2) >> 2;
            }
        }
    }

    /*Smooth an (height x width) image inplace with (121)>>2 vertically and horizontally.
    No access is performed outside the image, instead, we put special cases at the borders to avoid accessing outside line
    This function needs (2 x width) entries of scratch space in ss.*/
    inline void i_binomial_121_inplace(unsigned char **image, 
        int width, int height, 
        unsigned char **ss)
    {
        int i, im1, hm1;
        assert(height > 1 && width > 1);
        /*first row*/
        i_binomial_v_121(image[0], image[0], image[1], ss[0], width);
        
        hm1 = height - 1;
        for (i = 1; i < hm1; i++)
        {
            im1 = i - 1;
            i_binomial_v_121(image[im1], image[i], image[i + 1], ss[i & 1], width);
            i_binomial_h_121(ss[im1 & 1], image[im1], width);
        }
        
        /*second last row*/
        im1 = i - 1;
        i_binomial_v_121(image[im1], image[i], image[i], ss[i & 1], width);
        i_binomial_h_121(ss[im1 & 1], image[im1], width);
        
        /*last row*/
        i++;
        im1 = i - 1;
        i_binomial_h_121(ss[im1 & 1], image[im1], width);
    }

    /*Run (121)*0.25 on n entries of s and store result in n (n must be > 1) entries of d (which could be the same as s).
    Notice that s[-1] and s[n] will Not be accessed, instead, we put special cases at the borders to avoid accessing outside line.*/
    inline void i_binomial_h_121(const float *s, float *d, int n)
    {
        int i, count, nm;
        float a, b, c;
    
        if (n>1)
        {
            nm = n - 1;
            a = s[1];
            b = *s++;
            count = (nm >> 4);
            for (i = 0; i<count; i++)
            {
                c = *s++; 
                *d++ = (a + (b*2.0f) + c)*0.25f;
                a = *s++; 
                *d++ = (b + (c*2.0f) + a)*0.25f;
                b = *s++; 
                *d++ = (c + (a*2.0f) + b)*0.25f;
    
                c = *s++; 
                *d++ = (a + (b*2.0f) + c)*0.25f;
                a = *s++; 
                *d++ = (b + (c*2.0f) + a)*0.25f;
                b = *s++; 
                *d++ = (c + (a*2.0f) + b)*0.25f;
    
                c = *s++; 
                *d++ = (a + (b*2.0f) + c)*0.25f;
                a = *s++; 
                *d++ = (b + (c*2.0f) + a)*0.25f;
                b = *s++; 
                *d++ = (c + (a*2.0f) + b)*0.25f;
    
                c = *s++; 
                *d++ = (a + (b*2.0f) + c)*0.25f;
                a = *s++; 
                *d++ = (b + (c*2.0f) + a)*0.25f;
                b = *s++; 
                *d++ = (c + (a*2.0f) + b)*0.25f;
    
                c = *s++; 
                *d++ = (a + (b*2.0f) + c)*0.25f;
                a = *s++; 
                *d++ = (b + (c*2.0f) + a)*0.25f;
                b = *s++; 
                *d++ = (c + (a*2.0f) + b)*0.25f;
    
                c = *s++; 
                *d++ = (a + (b*2.0f) + c)*0.25f;
    
                a = b; 
                b = c;
            }
            count = nm & 15;
            for (i = 0; i < count; i++)
            {
                c = *s++; 
                *d++ = (a + (b*2.0f) + c)*0.25f;
                a = b; 
                b = c;
            }
            *d++ = (a + (b*2.0f) + a)*0.25f;
        }
    }

    /*Run (121)*0.25 to n entries in (s1,s2,s3) and store result in d*/
    inline void i_binomial_v_121(const float *s1, const float *s2, const float *s3, float *d, int n)
    {
        int i, c;
    
        c = (n >> 4);
        for (i = 0; i < c; i++)
        {
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
    
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
    
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
    
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
        }
        if (n & 15)
        {
            if (n & 8)
            {
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
    
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            }
            if (n & 4)
            {
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            }
            if (n & 2)
            {
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            }
            if (n & 1)
            {
                *d++ = ((*s1++) + ((*s2++)*2.0f) + (*s3++))*0.25f;
            }
        }
    }

    /*Smooth an (height x width) image inplace with (121)*0.25 vertically and horizontally.
    No access is performed outside the image, instead, we put special cases at the borders to avoid accessing outside line
    This function needs (2 x width) entries of scratch space in ss.*/
    inline void i_binomial_121_inplace(float **image, int width, int height, float **ss)
    {
        int i, im1, hm1;
        assert(height > 1 && width > 1);
        /*first row*/
        i_binomial_v_121(image[0], image[0], image[1], ss[0], width);
    
        hm1 = height - 1;
        for (i = 1; i < hm1; i++)
        {
            im1 = i - 1;
            i_binomial_v_121(image[im1], image[i], image[i + 1], ss[i & 1], width);
            i_binomial_h_121(ss[im1 & 1], image[im1], width);
        }
    
        /*second last row*/
        im1 = i - 1;
        i_binomial_v_121(image[im1], image[i], image[i], ss[i & 1], width);
        i_binomial_h_121(ss[im1 & 1], image[im1], width);
    
        /*last row*/
        i++;
        im1 = i - 1;
        i_binomial_h_121(ss[im1 & 1], image[im1], width);
    }

    /*Run (161)>>3 on n entries of s and store result in n (n must be > 1) entries of d (which could be the same as s).
    Notice that s[-1] and s[n] will Not be accessed, instead, we put special cases at the borders to avoid accessing outside line.*/
    inline void i_binomial_h_161(const unsigned char *s, unsigned char *d, int n)
    {
        int i, count, nm1 = n - 1;
        unsigned char a, b, c; /*tmp memorty to store intermidiate results*/
        assert(n > 1);
        a = s[1];
        b = *s++;
        count = (nm1 >> 4);
        for (i = 0; i < count; i++)
        {
            /*add 4 to handle roundoff error*/
            c = *s++; 
            *d++ = (a + (b << 1) + (b << 2) + c + 4) >> 3; 
            a = *s++; 
            *d++ = (b + (c << 1) + (c << 2) + a + 4) >> 3;
            b = *s++; 
            *d++ = (c + (a << 1) + (a << 2) + b + 4) >> 3;
    
            c = *s++; 
            *d++ = (a + (b << 1) + (b << 2) + c + 4) >> 3;
            a = *s++; 
            *d++ = (b + (c << 1) + (c << 2) + a + 4) >> 3;
            b = *s++; 
            *d++ = (c + (a << 1) + (a << 2) + b + 4) >> 3;
    
            c = *s++; 
            *d++ = (a + (b << 1) + (b << 2) + c + 4) >> 3;
            a = *s++; 
            *d++ = (b + (c << 1) + (c << 2) + a + 4) >> 3;
            b = *s++; 
            *d++ = (c + (a << 1) + (a << 2) + b + 4) >> 3;
    
            c = *s++; 
            *d++ = (a + (b << 1) + (b << 2) + c + 4) >> 3;
            a = *s++; 
            *d++ = (b + (c << 1) + (c << 2) + a + 4) >> 3;
            b = *s++; 
            *d++ = (c + (a << 1) + (a << 2) + b + 4) >> 3;
    
            c = *s++; 
            *d++ = (a + (b << 1) + (b << 2) + c + 4) >> 3;
            a = *s++; 
            *d++ = (b + (c << 1) + (c << 2) + a + 4) >> 3;
            b = *s++; 
            *d++ = (c + (a << 1) + (a << 2) + b + 4) >> 3;
    
            c = *s++; 
            *d++ = (a + (b << 1) + (b << 2) + c + 4) >> 3;
    
            a = b; 
            b = c;
        }
        count = nm1 & 15; /*handle left columns*/
        for (i = 0; i < count; i++)
        {
            c = *s++; 
            *d++ = (a + (b << 1) + (b << 2) + c + 4) >> 3;
            a = b; 
            b = c;
        }
        *d++ = (a + (b << 1) + (b << 2) + a + 4) >> 3;
    }

    /*Run (161)>>3 to n entries in (s1,s2,s3) and store result in d*/
    inline void i_binomial_v_161(const unsigned char *s1, 
        const unsigned char *s2, 
        const unsigned char *s3, 
        unsigned char *d, 
        int n)
    {
        int i, c;
        c = (n >> 4);
        for (i = 0; i<c; i++)
        {
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;

            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;

            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;

            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
            *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
            s2++;
        }
        if (n & 15)
        {
            if (n & 8)
            {
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
    
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
            }
            if (n & 4)
            {
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
            }
            if (n & 2)
            {
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
            }
            if (n & 1)
            {
                *d++ = ((*s1++) + ((*s2) << 2) + ((*s2) << 1) + (*s3++) + 4) >> 3;
                s2++;
            }
        }
    }

    /*Smooth an (height x width) image inplace with (161)>>3 vertically and horizontally.
    No access is performed outside the image, instead, we put special cases at the borders to avoid accessing outside line
    This function needs (2 x width) entries of scratch space in ss.*/
    inline void i_binomial_161_inplace(unsigned char **image, 
        int width, int height, 
        unsigned char **ss)
    {
        int i, im1, hm1;
        assert(height > 1 && width > 1);
        /*first row*/
        i_binomial_v_161(image[0], image[0], image[1], ss[0], width);
    
        hm1 = height - 1;
        for (i = 1; i < hm1; i++)
        {
            im1 = i - 1;
            i_binomial_v_161(image[im1], image[i], image[i + 1], ss[i & 1], width);
            i_binomial_h_161(ss[im1 & 1], image[im1], width);
        }
    
        /*second last row*/
        im1 = i - 1;
        i_binomial_v_161(image[im1], image[i], image[i], ss[i & 1], width);
        i_binomial_h_161(ss[im1 & 1], image[im1], width);
    
        /*last row*/
        i++;
        im1 = i - 1;
        i_binomial_h_161(ss[im1 & 1], image[im1], width);
    }

    /*Run (161)*0.125 on n entries of s and store result in n (n must be > 1) entries of d (which could be the same as s).
    Notice that s[-1] and s[n] will Not be accessed, instead, we put special cases at the borders to avoid accessing outside line.*/
    inline void i_binomial_h_161(const float *s, float *d, int n)
    {
        int i, count, nm;
        float a, b, c;
    
        if (n>1)
        {
            nm = n - 1;
            a = s[1];
            b = *s++;
            count = (nm >> 4);
            for (i = 0; i<count; i++)
            {
                c = *s++; 
                *d++ = (a + (b*6.0f) + c)*0.125f;
                a = *s++; 
                *d++ = (b + (c*6.0f) + a)*0.125f;
                b = *s++; 
                *d++ = (c + (a*6.0f) + b)*0.125f;
    
                c = *s++; 
                *d++ = (a + (b*6.0f) + c)*0.125f;
                a = *s++; 
                *d++ = (b + (c*6.0f) + a)*0.125f;
                b = *s++; 
                *d++ = (c + (a*6.0f) + b)*0.125f;
    
                c = *s++; 
                *d++ = (a + (b*6.0f) + c)*0.125f;
                a = *s++; 
                *d++ = (b + (c*6.0f) + a)*0.125f;
                b = *s++; 
                *d++ = (c + (a*6.0f) + b)*0.125f;
    
                c = *s++; 
                *d++ = (a + (b*6.0f) + c)*0.125f;
                a = *s++; 
                *d++ = (b + (c*6.0f) + a)*0.125f;
                b = *s++; 
                *d++ = (c + (a*6.0f) + b)*0.125f;
    
                c = *s++; 
                *d++ = (a + (b*6.0f) + c)*0.125f;
                a = *s++; 
                *d++ = (b + (c*6.0f) + a)*0.125f;
                b = *s++; 
                *d++ = (c + (a*6.0f) + b)*0.125f;
    
                c = *s++; 
                *d++ = (a + (b*6.0f) + c)*0.125f;
    
                a = b; 
                b = c;
            }
            count = nm & 15;
            for (i = 0; i<count; i++)
            {
                c = *s++; 
                *d++ = (a + (b*6.0f) + c)*0.125f;
                a = b; 
                b = c;
            }
            *d++ = (a + (b*6.0f) + a)*0.125f;
        }
    }

    /*Run (161)*0.125 to n entries in (s1,s2,s3) and store result in d*/
    inline void i_binomial_v_161(const float *s1, const float *s2, const float *s3, float *d, int n)
    {
        int i, c;
    
        c = (n >> 4);
        for (i = 0; i<c; i++)
        {
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
    
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
    
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
    
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
        }
        if (n & 15)
        {
            if (n & 8)
            {
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
    
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            }
            if (n & 4)
            {
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            }
            if (n & 2)
            {
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            }
            if (n & 1)
            {
                *d++ = ((*s1++) + ((*s2++)*6.0f) + (*s3++))*0.125f;
            }
        }
    }

    /*Smooth an (height x width) image inplace with (161)*0.125 vertically and horizontally.
    No access is performed outside the image, instead, we put special cases at the borders to avoid accessing outside line
    This function needs (2 x width) entries of scratch space in ss.*/
    inline void i_binomial_161_inplace(float **image, int width, int height, float **ss)
    {
        int i, im1, hm1;
        assert(height > 1 && width > 1);
        /*first row*/
        i_binomial_v_161(image[0], image[0], image[1], ss[0], width);
    
        hm1 = height - 1;
        for (i = 1; i < hm1; i++)
        {
            im1 = i - 1;
            i_binomial_v_161(image[im1], image[i], image[i + 1], ss[i & 1], width);
            i_binomial_h_161(ss[im1 & 1], image[im1], width);
        }
    
        /*second last row*/
        im1 = i - 1;
        i_binomial_v_161(image[im1], image[i], image[i], ss[i & 1], width);
        i_binomial_h_161(ss[im1 & 1], image[im1], width);
    
        /*last row*/
        i++;
        im1 = i - 1;
        i_binomial_h_161(ss[im1 & 1], image[im1], width);
    }

} /* namespace idl */

#endif